Phase diagram of harmonically confined one-dimensional fermions with attractive and 

repulsive interactions 
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We construct the complete U-^ phase diagram for harmonically confined ultracold fermionic 
atoms with repulsive and attractive interactions, (/i is the chemical potential and U the interaction 
strength.) Our approach is based on density-functional theory, and employs analytical expressions 
for the kinetic and correlation energy functionals. For repulsive interactions our calculations confirm 
previous numerical studies, and complement these by closed expressions for all phase boundaries and 
characteristic lines of the phase diagram, and by providing an explanation and solution for difficulties 
encountered in earlier density-functional work on the same system. For attractive interactions we 
propose a new and accurate interpolation for the correlation energy, and use it to extend the phase 
diagram to ?7 < 0. 
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The field of strongly correlated low-dimensional 
fermions has a long history and multifarious applica- 
tions in condensed-matter physics. More recent exper- 
imental progress in the field of confined ultracold atoms 
Q, 1^, 1^ |j{ has further widened the scope of this already 
very rich field by introducing, e.g., an unprecedented de- 
gree of control over the confining potential, and the pos- 
sibility to continuously tune the particle-particle interac- 
tion from repulsive to attractive — which is not easy for 
cold atoms, but unconceivable for electrons. After much 
progress on bosonic atoms and their condensation, more 
recently fermionic atoms have come under intense study. 

In the present paper we focus on one particular 
model for optically confined fermionic atoms, the one- 
dimensional Hubbard model. The parameters character- 
izing this model and spanning its phase diagram are the 
on-site interaction U ^ and the filling factor n = N/L, 
where N is the number of particles and L the number 
of lattice sites. For systems in which n is spatially vary- 
ing, such as in the presence of a confining potential, it is 
more useful to specify the chemical potential /i instead, 
which is a constant throughout the system. Hence, we are 
interested in the phase diagram of the one-dimensional 
Hubbard model in U-fi space. Our aim is to map out 
this phase diagram, identify the physics governing each 
possible phase, and to provide analytical expressions for 
all phase boundaries. 

The U > 0, i.e., repulsive, part of this phase diagram 
has recently been studied also by Liu et al. 0, who used 
a half-numerical half-analytical local-density approxima- 
tion (LDA). Rigol et al. [3 obtained essentially the same 
U > phase diagram from fully numerical Quantum 
Monte Carlo simulations. An interesting feature that 
emerged in these works is the appearance of mixed Mott- 
insulator Luttinger-liquid phases, characterized by spa- 



tially separated incompressible and compressible regions 
in the optical trap, and plateaus in the density profiles 
coinciding with the incompressible regions. Xianlong et 
al. have recently studied the case of negative U, i.e., 
attractive interactions, where the physics is dominated by 
Luther-Emery-type pairing correlations and the sponta- 
neous formation of atomic density waves. However, the 
Kohn-Sham-type density-functional calculations of Ref. 
do not allow one to map out the entire phase diagram, 
because the self-consistency cycle does not converge in 
the mixed metal-insulator phases aX U > 0. 

In the present paper we use an analytical local-density 
approximation both for the kinetic energy (in the spirit of 
Thomas-Fermi theory) and for the correlation energy (in 
the spirit of the Bethe-Ansatz LDA proposed in Ref. |9j) 
to obtain the U-n phase diagram. For repulsive inter- 
actions our results coincide with those obtained numeri- 
cally in Refs. HQ. However, as our approach is entirely 
analytical, we obtain closed expressions for the various 
phase boundaries, which previously had to be extracted 
numerically. Moreover, our present work also covers the 
case of attractive interactions, t7 < 0, thus mapping out 
the entire physically relevant phase diagram. 

To make our paper self-contained we start by briefly 
describing and comparing the various approximation 
schemes employed. Density-functional theory (DFT) 
[Tol ITU shows that the total energy of a many-body sys- 
tem in an arbitrary potential can be written 



E[n] = Ts[n] + EhM + V[n] + E.,,[r 



(1) 
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where Ts is the kinetic energy of noninteracting particles, 
Eh is the Hartree energy (i.e., the mean-field approx- 
imation to the interaction energy), V is the potential 
energy (arising from the spatial distribution of the nuclei 
in solid-state and molecular applications of DFT, and 
from the optical confining potential for atoms in optical 
traps), and E^c is the exchange-correlation energy _12|. 
The local-density approximation to any energy compo- 



nent F is defined by F 
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tinuum, and F^'^^ = X^i on a lattice, where 
/ denotes the per-volume or per-site value of F in a spa- 
tially uniform (homogeneous) system. 

Although DFT is mostly applied to ab initio electronic- 
structure calculations in solid-state physics and quantum 
chemistry [ToL [Tll |. its formal framework is sufficiently 
general to permit application also to wide classes of 
model Hamiltonia ns »13il . and, in particular, to the Hub- 
bard model IEQIIJa Bethe-Ansatz based LDA func- 
tional for the correlation energy of the one-dimensional 
Hubbard model was proposed in and applied to Mott 
insulating phases 0| and superlattice structures \it\ . 
These works employed a parametrization of the per-site 
ground-state energy of the spatially uniform Hubbard 
model, e(n, U), which was constructed to recover ex- 
actly known results at J7 = 0, C/ — > oo and n = 1, and to 
be close to data obtained from numerical solution of the 
Bethe-Ansatz (BA) equations inbetween jl^. For 
< n < 1 



2/3(t/) . f 7:n \ 



and for 1 < n < 2 



<n,U) = {r 



2(3(U) . (ir{2 
l)U sm 



(2) 
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Here I3{U) is obtained from 



/3 . 



= -2 



Jo{x)Ji{x) 
a;(l + e^^^/2) 



dx 



-2/(f/), 



(4) 

and Jm is the m'th order Bessel function. The difference 
between the n < 1 and the n > 1 branch is responsible 
for a derivative discontinuity of Ec and the opening of a 
Mott gap at n = 1 [H 18]. 

To perform analytical calculations for attractive inter- 
actions, our first task is to extend this construction to 
U < 0. The construction of a parametrization of the per- 
site energy is not unique, but a convenient and accurate 
choice is 

e(n,t/) = ^-4/(|C/|)sin(^). (5) 

This expression by construction satisfies the following ex- 
actly known properties of e(n, U) for negative U: (i) 
e(n = 0,[/) = 0, (ii) e{n ^ l,U) ^ ~4/(|C/|) -I- U/2, 
(iii) e{n,U = 0) = -(4/7r) sin(7rn/2), (iv) e{n,U 
-oo) ^ Un/2, (v) de{n,U)/dn\n=i = U/2, and (vi) 
e(2 -n,U) + U{n - 1) = e(n, U). Note that unlike @ 
and (|2Jl, the U < expression (jSJ is particle- hole sym- 
metric [property (iv)], and defined on only one branch. 

The Kohn-Sham approach treats Ts exactly, by means 
of single-particle orbitals, and approximates Ec, e.g., by 
the BA-LDA of Ref. M. 



Eks U] = Ts[n] -f V[n] + Enin, U] 



where for a spin-unpolarized system EH[n,U] = 
(C//4)X- nf, and V[n] = J2i''^i^i the potential en- 
ergy in the confining potential Vi. ts{n) — e(n,U = 0) 
is the per-site kinetic energy of noninteracting particles, 
and the difference e(n, U) — e//(n, U) — ts{n) =: edn, U) 
is by definition the correlation energy. 

The total-energy LDA (TLDA), employed below, uses 
an LDA both for the kinetic and the correlation energy, 
so that 

iJ^^^^[n] =^t,(n,) + yW+i?HW+^ec(n,) (7) 

i i 

= ^e(n„C/) + F[n], (8) 

i 

where the second equality follows because the Hartree en- 
ergy is a local functional for an on-site interaction. Note 
that the TLDA approach is different from the Thomas- 
Fermi approximation (TFA), which also employs an LDA 
for Ts, but takes Ec = 0. 

We now proceed to characterize the phase diagram. To 
this end we employ the TLDA, which is more accurate 
than the TFA, and unlike KS procedures can yield closed 
analytical results. The density profile is obtained from 
minimizing the energy functional under the constraint of 
fixed total particle number. Hence 



d{E-nN) _ de{ni,U) 
Srij duj 



(9) 



From the above parametrizations © and jSJ of 
e(n, U) we obtain the Euler equations 



2 cos ] + fi = < 77 i < 1 (10) 

[/4-2cos( ^ " \ + V - fi^O l<ni<2(ll) 



for U > 0, and 
U 



--27r/(lC/l)cos(^)+l/.-M = ( 



12) 



for U < 0. Note that at = 1 the expression e{n, U) 
for [/ > is not differentiable and no Euler equation is 
obtained. The three Euler equations (|10() - (|12l) have a 
complex solution space with a rich structure, which we 
now proceed to analyse. 

Let us first consider f/ > 0. In the low-density region 
far from the trap center the solution of Eq. (|10|l is rii — 

^ arccos ^ ^'2^^ ^ . However, this solution only belongs to 

the interval [0, 1] if 

V - 2 < fi <V - 2cos (n/P) . (13) 

Closer to the center of the trap Ui > 1, and the rel- 
evant branch of the solution is Hll|) . leading to Ui — 
2 — (P/tt) arccos((/^ — Vi — U)/2), provided that 

V + U + 2cos{7t/(3) < fi<V + U + 2. (14) 
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FIG. 1: Classification of density sliapes across tfie i > lialf 
of a harmonic trap with Vi = 0.004i^, for [/ = 4 > (7*. The 
four dividing Unes have the form Vi + A, where, from top to 
bottom, A = U + 2, A = U + 2co8{n/ f3), A = -2cos(7r//3) 
and A = —2. The horizontal line at = 4 is the example 
discussed in the main text. 

Note that the conditions and (|14|) . which guarantee 
that the solution found does indeed pertain to the den- 
sity interval it is obtained from, do not necessarily match 
continuously. In particular, if 

C/ + 2 cos (tt//?) > -2 cos {n/(3) (15) 

there are values of fi for which neither of the two branches 
of the Euler equation has a solution. The actual density 
in this situation is then the one for which the Euler equa- 
tions are not defined, i.e., Ui = 1. This possibility occurs 
only for U > U* , where U* = 1.7349... is the lowest value 
of U satisfying ((TB|l . 

In Fig. n we plot the four limiting expressions defining 
the boundaries of Eqs. l(T3Jl and lO for U > U*. The 
chemical potential is constant throughout the trap, and 
thus corresponds to horizontal lines crossing the various 
boundaries. As an illustration, consider a system with 
/i = 4 (horizontal line in Fig. . At the trap center {i = 
0) one first finds a high-density solution with 1 < rii < 2. 
Moving outward, one crosses the left boundary of (|14|l . 
encountering a region in which neither Euler equation has 
a solution and the density is fixed at = 1 . Further out 
a solution of H13|l becomes possible, and for z > 39 even 
the low-density Euler equation ceases to have a solution, 
implying Ui = 0. Density profiles for [7 = 4 and different 
fillings are shown in Fig. [3 Each of the data sets in Fig. [21 
illustrates one of the five phases of the U > part of the 
phase diagram, Fig. |21 These are the same five phases 
also obtained numerically in Refs. 0, Q . 

The present TLDA approach localizes the sites with 
rii = 1 by the conditions that < n < 2 and fj, such that 
neither Euler equation has a solution. On the other hand, 
m a Kohn-Sham approach H Hill [13 a self-consistent 
solution is obtained iteratively, yielding in each iteration 
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FIG. 2: Typical density profiles for (7 = 4 and a confining 
potential Vi = 0.004j^, for five different values of the total 
fermion number, ilustrating the five basic types of density 
profiles obtained for repulsive interactions. 

step a density that is defined simultaneously at all sites. 
If for particular sites there is no value of rii corresponding 
to a solution of the Euler equations, the Kohn-Sham self- 
consistency cycle does not find any solution at all and 
does not converge. This is precisely what was observed 
in trying to perform calculations of the type of Refs. [1,13 
in a region with plateaus. The TLDA procedure, which 
allows to obtain and characterize solutions site by site, 
has no such problem. 

The negative U analysis is performed in the same way, 
by starting from Eq. H12() . Since there is only one branch, 
the criterium for a valid solution is 

^ + V^- 2nIi\U\) <^i<J + V, + 2TrI{\U\). (16) 

For n inside this interval the solution of H12|l is rii = 
{2/^T)a.rccos{{V^ + U/2- fi)/{2TrI{\U\))). For ^ < f -|- 
Vi ~ 2ttI{\U\) the solution is = and for fi > ^ + 
Vi + 27r/(|C7|) it is rii — 2. Typical density profiles in this 
region (not shown) are of the same type as for [7 > 0, 
with possible plateaus at n = and n — 2. The absence 
of any possible plateaus with n = 1 is due to the absence 
of a derivative discontinuity of Ec. The atomic density 
waves found from Kohn-Sham calculations in Ref. [H are 
not reproduced by the TLDA, showing that observation 
of these oscillations (just as Friedel oscillations) require 
an exact (orbital) treatment of the single-body kinetic 
energy. 

The expresions for the phase boundaries obtained 
within TLDA can be used to analytically construct the 
full U-^ phase diagram, shown in Fig. 13 For a harmonic 
confining potential of the form Vi — ki^ region II (char- 
acterized by a central plateau with n — 1, i.e., a local 
Mott-insulator-like state) is bounded by fc — 2 cos(7r//3) < 
^1 < U + 2cos(7r//3). Region III (two lateral plateaus 
with n = 1) is bounded by f7 -|- 2cos(7r//3) < ^ < 
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FIG. 3: The (7-/1 phase diagram. See main text for details. 

k + U/2 + [U/A + coa{n/l3))^/k, and region V (a central 
band-insulator-like plateau with n = 2) by > U + 2 + k. 
Region IV is the set of points satisfying simultaneously 
the criteria for region III and V, i.e., has lateral plateaus 
with n — 1 and a central plateau with n = 2. Region I 
(purely metallic without any plateaus, except for, possi- 
bly, sites with n = at the wings of the trap) is the set 
of points with /i > — 2 that do not belong to any other 
U > region. In the TFA {Ec = 0) the f7 > phases 
II, III and IV disappear, and the dividing line between 
phase I and V becomes fi = U + 2 + k. Phases II, III and 
IV thus owe their existence to correlation. Region VI is 
the U < extension of region V and characterized by a 
plateau at n = 2 in the trap center. It is bounded from 
below hy k + U/2 + 27r/(|J7|). Region VII, finally, is 



the 17 < extension of region I and characterized by the 
absence of any plateaus except for, possibly, sites with 
n = at the wings of the trap. 

Special values of U and /i can be used to further classify 
the possibilities: For < /i*, where fi*{U > 0) = — 2 
and IJ.*{U < 0) = U/2 - 2ttI{\U\), the only solution is 
rii = at all sites, i.e., an empty system (hatched region 
in the phase diagram). For < U < U* the intervals 
permitting solutions of H13|l and H14|) overlap, and no n = 
1 plateaus can appear, leaving only phases I, V and VI. 
(In this situation, at isolated sites, an rt > 1 and an 
n < 1 solution can coexist.) For U* < U < U**{k), where 
U**{k) is the solution of {U/4, + cos{-K/l3) f = k{U/2 + 2), 
there cannot be a phase IV (i.e, plateaus at n = 1 and 
n — 2) and instead there appears a reentrant metallic 
phase of type I above phase III. 

In summary, we have constructed the complete U-^ 
phase diagram of one-dimensional interacting harmon- 
ically confined fermions. For [/ > we confirm earlier 
numerical results, but complement them by providing an- 
alytical expressions for the phase boundaries. The TLDA 
approach also provides an explanation why Kohn-Sham 
calculations do not work in the metal Mott-insulator 
phase-separated state (with plateaus at n = 1) and pro- 
vides an alternative way for obtaining density profiles in 
this region. For J7 < we propose, in Eq. (Q, a simple 
analytical parametrization of the Bethe-Ansatz solution, 
which we use, within TLDA, to extend the phase diagram 
to attractive interactions. 
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